Optimization by Quantum Annealing: Lessons from Simple Cases 
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This paper investigates the basic behavior and performance of simulated quantum annealing (QA) 
in comparison with classical annealing (CA). Three simple one dimensional case study systems are 
considered, namely a parabolic well, a double well, and a curved washboard. The time dependent 
Schrodinger evolution in either real or imaginary time describing QA is contrasted with the Fokker 
Planck evolution of CA. The asymptotic decrease of excess energy with annealing time is studied 
in each case, and the reasons for differences are examined and discussed. The Huse-Fisher classical 
power law of double well CA is replaced with a different power law in QA. The multi-well wash- 
board problem studied in CA by Shinomoto and Kabashima and leading classically to a logarithmic 
annealing even in the absence of disorder, turns to a power law behavior when annealed with QA. 
The crucial role of disorder and localization is briefly discussed. 
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I. INTRODUCTION 

The idea of quantum annealing (QA) is a late off- 
spring of the celebrated simulated thermal annealing by 
Kirkpatrick et aim. In simulated annealing, the prob- 
lem of minimizing a certain cost (or energy) function in 
a large configuration space is tackled by the introduc- 
tion of a fictitious temperature, which is slowly lowered 
in the course of a Monte Carlo or Molecular Dynamics 
simulation^. This device allows an exploration of the 
configuration space of the problem at hand, effectively 
avoiding trapping at unfavorable local minima through 
thermal hopping above energy barriers. It makes for a 
very robust and effective minimization tool, often much 
more effective than standard, gradient-based, minimiza- 
tion methods. 

An elegant and fascinating alternative to such a clas- 
sical simulated annealing (CA) consists in helping the 
system escape the local minima through quantum me- 
chanics, by tunneling through the barriers rather than 
thermally overcoming thern^. Experimental evidence in 
disordered Ising ferromagnets subject to transverse mag- 
netic fields showed that this strategy is not only feasible 
but presumably winning in certain casesi In essence, 
in quantum annealing one supplements the classical en- 
ergy function, let us denote it by H c i, with a suitable 
time- dependent quantum kinetic term, Hkm(t), which is 
initially very large, for t < 0, then gradually reduced to 
zero in a time r. The quantum state of the system |\I/(t)), 
initially prepared in the fully quantum ground state \^o) 
of H(t = 0) = H c i + i?fei„(0), evolves according to the 
Scrodinger equation 



m-Mt)} 



[H d + H kin {t)]\V(t)) , 



(1) 



to reach a final state \^(t = t)). A crucial basic question 
is then how the residual energy e res (r) = Efi n {r) — E opt , 
decreases for increasing r. Here E opt is the absolute min- 



imum of H c i, and Efi n (r) is the average energy attained 
by the system after evolving for a time r, Efi n (r) — 
(*(r)|if c {|*(r)>/(*(r)|*(r)>, Generally speaking, this 
question has to do with the adiabaticity of the quantum 
evolution, i.e., whether the system is able, for sufficiently 
slow annealing (sufficiently long r), to follow the instan- 
taneous ground state of H(t) = H c i + Hun{t), for a ju- 
diciously chosen H^ n {t). The fictitious kinetic energy 
Hkin (t) can be chosen quite freely, with the only require- 
ment of being reasonably easy to implement. For this 
reason, this approach has also been called Quantum Adi- 
abatic Evolution^. 

At the level of practical implementations on an ordi- 
nary (classical) computer, the task of following the time- 
dependent Schrodinger evolution in Eq. ^ is clearly fea- 
sible only for toy models with a sufficiently manageable 
Hilbert spaced. Actual optimization problems of practi- 
cal interest usually involve astronomically large Hilbert 
spaces, a fact that calls for alternative stochastic (Quan- 
tum Monte Carlo) approaches. These QMC techniques, 
in turn, are usually suitable to using imaginary time 
quantum evolution, where the ihdt in Eq. ^ is replaced 
by —hdt- One of the questions we will try to address in 
the present paper will be, therefore, if working with an 
imaginary-time Schrodinger evolution changes the quan- 
tum adiabatic evolution approach in any essential way. 
We will argue that, as far as annealing is concerned, 
imaginary-time is essentially equivalent to real-time, and, 
as a matter of fact, can be quantitatively better. 

Alternatively, a number of recent theoretical papers 
have applied Path-Integral Monte Carlo strategies to QA. 
A certain success has been obtained in several optimiza- 
tion problems, such as the folding of off- lattice poly- 
mer models^* 7 , the random Ising model ground state 
problem^, and the Traveling Salesman ProblemiS. It 
is fair to say, however, that there is no general theory 
predicting the performance of a QA algorithm, in par- 
ticular correlating it with the energy landscape of the 
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given optimization problem. This is a quite unpleasant 
situation, in view of the fact that it is a-priori not obvi- 
ous or guaranteed that a QA approach should do better 
than, for instance, CA. Indeed, for the interesting case of 
Boolean Satisfiability problems - more precisely, a pro- 
totypical NP-complete problem such as 3-SAT - recent 
attempts in our group showed that Path-Integral Monte 
Carlo annealing may perform definitely worse than sim- 
ple CAii. 

Evidently, the performance of QA over CA depends 
in detail on the energy landscape of the problem at 
hand, in particular on the nature and type of barriers 
separating the different local minima, a problem about 
which very little is known in many practical interesting 
casesi^. It makes sense therefore to move one step back 
and concentrate attention on the simplest textbook prob- 
lems where the energy landscape is well under control: 
essentially, one-dimensional potentials, starting from a 
double well potential, the simplest form of barrier. On 
these well controlled landscapes we can carry out a de- 
tailed and exhaustive comparison between quantum adi- 
abatic Schrodinger evolution, both in real and in imag- 
inary time, and its classical deterministic counterpart, 
i.e., Fokker-Planck evolution. 

The rest of the paper is organized as follows. In Sec.llTI 
we define the problem we want to tackle, i.e., comparing 
Fokker-Planck annealing to Schrodinger annealing. In 
Sec. IIIII we consider in detail the case of a double-well 
barrier. In Sec. II VI we move to the case of a potential 
with many minima, but no disorder, were the behavior of 
classical and quantum annealing is remarkably different. 
In Sec.|V]the crucial role of disorder is discussed. Finally, 
Sec. IVII contains a summary of the results found and a 
few concluding remarks. Details of the calculations are 
contained in the Appendices. 



II. SCHRODINGER VERSUS 
FOKKER-PLANCK ANNEALING: STATEMENT 
OF THE PROBLEM 

Suppose we are given a potential V(x), (with x a Carte- 
sian vector of arbitrary dimension) , of which we need to 
determine the absolute minimum (x op t,E opt = V(x op t)). 
Assume generally a situation in which a steepest-descent 
approach, i.e., the strategy of following the gradient of 
V, would lead to trapping into one of the many local 
minima of V , and would thus not work. Classically, as 
an obvious generalization of a steepest-descent approach, 
one could imagine of performing a stochastic (Markov) 
dynamics in x-space according to a Langevin equation 



(2) 



where the strength o f the noise term £ is controlled by the 
squared correlations = 2D(T)SijS{t — t'), with 

£ = 0. Both D(T) and T](T) - with dimensions of a diffu- 
sion constant and of a friction coefficient and related, re- 



spectively, to fluctuations and dissipation in the system - 
are temperature dependent quantities which can be cho- 
sen, for the present optimization purpose, with a certain 
freedom. The only obvious constraint is in fact that the 
correct thermodynamical averages should be recovered 
from the Langevin dynamics only if rj(T)D{T) = ksT, 
an equality known as Einstein's relation 13 . Physically, 
D(T) should be an increasing function of T, so as to 
lead to increasing random forces as T increases, with 
D(T = 0) = 0, since noise is turned off at T = 0. Clas- 
sical annealing can in principle be performed through 
this Langevin dynamics, by slowly decreasing the tem- 
perature T(t) as a function of time, from some initially 
large value To down to zero. Instead of working with the 
Langevin equation - a stochastic differential equation - 
one might equivalently address the problem by studying 
the probability density P(x, t) of finding a particle at po- 
sition x at time t. The probability density is well known 
to obey a deterministic time-evolution equation given by 
the Fokker-Planck (FP) equation^: 
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tj(T) 



div (PVV) + D(T)W 2 P 



(3) 



Here, the second term in the right-hand side represents 
the well known diffusion term, proportional to the diffu- 
sion coefficient D(T), whereas the first term represents 
the effect of the drift force — VV, inversely proportional 
to the friction coefficient n(T) = k B T/D(T)ll. Anneal- 
ing can now be performed, in principle, by keeping the 
system for a long enough equilibration time at a large 
temperature To, and then gradually decreasing T to zero 
as a function of time, T(t), in a given annealing time r. 
We can model this by assuming 



T(t) = T f(t/r) 



(4) 



where f(y) is some assigned monotonically decreasing 
function for y G [0, 1], with f(y < 0) = 1 and /(I) = 0. In 
this manner the diffusion constant D in Eq. JSJ becomes 
a time-dependent quantity, D t = D(T(t)). The FP equa- 
tion should then be solved with an initial condition given 
by the equilibrium Boltzmann distribution at tempera- 
ture T(t = 0) = T , i.e., P(x,t = 0) = e - v ^y kBTo . The 
final average potential energy after annealing in excess of 
the true minimum value, will then be simply given by: 

eres(T) = J dxV{x)P(x,t = r) - E opt > , (5) 

where E opt is the actual absolute minimum of the poten- 
tial V. 

In a completely analogous manner, we can conceive 
using Schrodinger 's equation to perform a determinis- 
tic quantum annealing (QA) evolution of the system, by 
studying: 

ih^l>{x, t) = [-T{t)V 2 + V{x)\ f/,(x, t) , (6) 

where £ = i for a real-time (RT) evolution, while £ = — I 
for an imaginary-time (IT) evolution. Here T(t) — 
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dt 



'pot 



(t) = kD t 



k B T(t) 



'-pot 



(*) 



h 2 /2mt will be our annealing parameter, playing the potential energy e po t(t), which has the form: 
role that the temperature T(t) had in classical anneal- 
ing. Once again we may take T(t) varying from some 
large value Tq sd t < - corresponding to a small mass 
of the particle, hence to large quantum fluctuations - 
down to T(t = r) = 0, corresponding to a particle of in- 
finite mass, hence without quantum fluctuations. Again, 
we can model this with T(t) — Tof(t/r), where / is a 
preassigned monotonically decreasing function. A conve- 
nient initial condition here will be tp(x,t = 0) = ipo( x )> 
where tpo(x) is the ground state of the system at t < 0, 
corresponding to the large value T(t) = T and hence 
to large quantum fluctuations. For such a large T the 
ground state will be isolated, and separated by a large 
energy gap from all excited states. The residual energy 
after annealing will be similarly given by Eq. (JSJ, where 
now, however, the probability P(x, t = r) should be in- 
terpreted as 



(7) 



the initial condition being simply given by the equipar- 
tition value e pot (t = 0) = ksT /2. As with every one- 
dimensional linear differential equation, Eq. (JJJ can be 
solved by quadrature for any choice of T(t) and D t = 
D(T(t)). Assuming the annealing schedule to be param- 
eterized by an exponent ar > 0, T(t) = Tb(l — i/r) QT , 
r being the annealing time, and the diffusion constant 
D(T) to behave as a power law of temperature, D(T) = 
Do(T/To) aD with «c > 0, we can easily extract from the 
analytical solution for e po t{t) the large r asymptotic be- 
havior of the final residual energy e res (r) = e pot (t = r) . 
That turns out to be: 



.-Sir 



with Qca — 



ar(aD - 1) + 1 



P(x,t) = 



|#M)| S 



/ dx' IV'O',*)! 5 



In general, the residual energy will be different for a RT 
or an IT Schrodinger evolution. We will discuss in some 
detail RT versus IT Schrodinger evolution in Sec. IIII Al 
The basic question we pose is which annealing scheme 
is eventually more effective, leading to the smallest final 
residual energies e res (r). This might seem an ill-posed 
question, because the time scales involved in the classi- 
cal and in the quantum evolution are different, and also 
because, practically, the two approaches might involve 
different computational costs which would imply differ- 
ent CPU time scales. In other words, it might seem that 
it only makes sense to ask how well an annealing scheme 
performs in a given CPU-time Tcpu, with all the un- 
avoidable uncertainty associated to a CPU-time-related 
answer (involving, among other things, the programmer's 
skills, the algorithmic choices, and the computer archi- 
tecture). We will show, however, that the behavior of 
£res( T ) can be so vastly different for the different schemes, 
obviously in strict relation with the form of the poten- 
tial, that such time scale concerns are often practically 
irrelevant 15 . 



A. The harmonic potential: a warm-up exercise 

Preliminary to any further treatment of a potential 
with barriers, and as a warm-up exercise which will be 
useful later on, we start here with the simple case of 
a parabolic potential in one dimension, V(x) — kx 2 /2, 
which has a trivial minimum in x = 0, with E opt = 0, 
and no barriers whatsoever. 

Let us consider classical FP annealing first. As de- 
tailed in the Appendix, it is a matter of simple algebra 
to show that, for the harmonic potential, one can write a 
simple closed linear differential equation^ for the average 



( 8 ) 

Trivial as it is, annealing proceeds here extremely fast, 
with a power-law exponent Qca that can increase with- 
out bounds (for instance if qb = 1) upon increasing the 
exponent cut of the annealing schedule. 

Consider now the Schrodinger evolution problem for 
this potential, 



d 



tp(x,t = 0) = ip {x) , 



i>(x,t) (9) 



where ipo(x) cx exp (—Bqx 2 /2) is the ground state Gaus- 
sian wavefunction corresponding to the initial value of 
the Laplacian coefficient T(t = 0) = r , and £ = ih or 
£ = — h for a real time (RT) or an imaginary time (IT) 
evolution, respectively. This problem is studied in detail 
in the Appendix, where we show that a Gaussian Ansatz 
for tp(x,t), of the form ip(x,t) oc exp (— Btx 2 /2) with 
Real(£?t) > 0, satisfies the time-dependent Schrodinger 
equation as long as the inverse variance B t of the Gaus- 
sian satisfies the following ordinary non-linear first-order 
differential equation: 



-£Bt - k-2T{t)B, 
B t= o = Bq = 



(10) 



k 
2T~ 



Contrary to the classical case, there is no simple way of 
recasting the annealing problem in terms of a closed lin- 
ear differential equation for the average potential energy 
£ P ot(t). The final residual energy e res (r) = e po t(t = t) is 
still expressed in terms of B T (or better, of its real part 
K(B T )), 



^res ("7") 



JdxV{x)\ip(x,t = T)\^ 
Jdx\ip{x,t = T)\ 2 



45R(B T ) 



(11) 



but the behavior of B t must be extracted from the study 
of the non- linear equation (|10|l . The properties of the so- 
lutions of Eq. I|1(J|) are studied in detail in the Appendix, 
where we show that: 
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i) eresfV) cannot decrease faster than 1/r, for large 
t, i.e., a power-law exponent e res (r) « r~ fi< ? A is 
bound to be < 1. 

ii) Adopting a power-law annealing schedule T(t) = 
Tq(1 — t/r) ar , the exponent Qqa for the IT case is 



the splitting between the two minima is given by Ay 



a r + 2 ' 



(12) 



increasing towards the upper bound 1 as ar is in- 
creased towards oo. 

iii) RT quantum annealing proceeds with exactly the 
same exponent Qqa as IT quantum annealing - al- 
though e^J(r) > e^(r) in general -, except that 
the limit ar — > oo (abrupt switch-off of the Lapla- 
cian coefficient) is singular in the RT case. 

Summarizing, we have learned that, for a single 
parabolic valley in configuration space, both CA and QA 
proceed with power-laws, but CA can be much more ef- 
ficient than QA, with an arbitrarily larger power-law ex- 
ponent. We underline however that this is merely an aca- 
demic matter at this point, steepest descent being much 
more efficient than both CA and QA in such a simple 
case. The power of QA shows up only when potentials 
with barriers are considered. 



III. THE SIMPLEST BARRIER: A 
DOUBLE- WELL POTENTIAL 

Let us take now the classical potential to be optimized 
as a simple double-well potential in one-dimension 



V sym (x) = V a 



(x 2 - a 2 ) 



2\2 



5x , 



(13) 



with Vo, a, and 5 real constants. In absence of the linear 
term (5 = 0), the potential has two degenerate minima 
located at ±a, and separated by a barrier of height Vo- 
When a small linear term 6 > is introduced , with 
5a << Vo, the two degenerate minima are split by a 
quantity Ay rs 25a, the minimum at x ~ —a becom- 
ing slightly favored. For reasons that will be clear in a 
moment, it is useful to slightly generalize the previous po- 
tential to a less symmetric situation, where the two wells 
possess definitely distinct curvatures at the minimum (i.e, 
their widths differ substantially). This is realized easily, 
with a potential of the form: 



Vasym (x) 



Vo 



V 



t 2 2 \2 

(x — a_) 



5x 
5x 



for x > 
for x < 



(14) 



with a+ 7^ a_, both positive. (The discontinuity in the 
second derivative at the origin is of no consequence in our 
discussion.) To linear order in the small parameter 5, the 
two minima are now located at x± = ±a± — Sa±/ (8Vq), 



5(a A 



.), while the second derivative of the potential 



at the two minima, to lowest order in 5, is given by: 

V (x = x±) = -5- . 
Obviously, V sym is recovered if we set a+ = a_ = a in 

Vasym ■ 

We now present the results obtained by the anneal- 
ing schemes introduced in Sec. ^ above. The Fokker- 
Planck and the Schrodinger equation (both in RT and 
in IT) were integrated numerically using a fourth-order 
Runge-Kutta method, after discretizing the x variable 
in a sufficiently fine real space gridit For the FP clas- 
sical annealing, the results shown are obtained with a 
linear temperature schedule, T(t) = To(l — t/r) (i.e., 
olt = 1), and a diffusion coefficient simply proportional 
to T{t), D t = D (l - t/r) (i.e., a D = 1). Conse- 
quently, the friction coefficient is kept constant in t, rj t = 
kBT(t)/D t = kBTo/Do- Similarly, for the Schrodinger 
quantum annealing we show results obtained with a co- 
efficient of the Laplacian T(t) vanishing linearly in a time 
t, T(t) = r (l - t/r) (i.e., a r = 1). 

Fig.^shows the results obtained for the final annealed 
probability distribution P(x, t = t) at different values 
of r, for both the Fokker-Planck (CA, panel (a)) and 
the Scrodinger imaginary-time case (IT, panel (b)), for 
a symmetric double well potential Ksym^), with Vo = 1 
(our unit of energy), a = a + = a_ = 1 (unit of length), 
5 = 0.1. Fig. ^c) summarizes the results obtained for 
the residual energy e res (T). Fig.|2Ia,b,c) shows the corre- 
sponding results for an asymmetric double well potential, 
Eq. HI with a+ = 1.25, a_ = 0.75, and 6 = 0.1. 

We notice immediately that QA wins, in both cases, 
over CA for large enough value of r. The RT-QA behaves 
as its IT counterpart for the symmetric double well, while 
it shows a different behavior in the asymmetric case (see 
below for comments). To go deeper into the details of 
the different evolutions, let us begin discussing the CA 
data (panel (a) and (c) of Figs. and EJ, which show 
similar behaviors for both choices of the potential. Start- 
ing from an initially broad Boltzmann distribution at a 
high T = T = Vo, P(x,t = 0) (solid lines), the system 
quickly sharpens the distribution P(x, t) into two well- 
defined and quite narrow peaks located around the two 
minima x± of the potential. This agrees very well with 
expectations based on the CA in a harmonic potential, 
which showed that the width of the Gaussian should de- 
crease linearly in r (flcA = I for ar = aD = I), as 
is indeed found in our double well case too. If we de- 
note by p± the integral of each of the two narrow peaks, 
with p_ = 1, it is clear that the problem has effec- 
tively been reduced to a discrete two-level system prob- 
lem. The time evolution p± therefore obeys a discrete 
Master equation which involves the thermal promotion 
of particles over the barrier Vo, of the form: 
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-i dp+ 
dt 



[l-p + (t)]e ^Srm- p+ (t)e : (15) 
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FIG. 1: (a,b): The annealed final probability distribution 
P(x,t = r) at different values of the annealing time r, for 
both the Fokker-Planck classical annealing (CA, panel (a)), 
and the Imaginary Time Schrodinger quantum annealing (IT- 
QA, panel (b)). (c) Final residual energy e re s(i~) versus an- 
nealing time r for quantum annealing in Real Time (RT) and 
Imaginary Time (IT) compared to the Fokker-Planck classical 
annealing (CA). The solid line in (c) is a fit of the CA data. 
The double well potential (dashed line in (a,b), inset of (c)) 
is here given by Eq. 1141 with a+ = a~ = a = 1. 



where 7 is an attempt frequency, while B = Vo — V{x+) 
and B + Ay = Vo — V(x- ) are the potential barriers seen 
from a particle in the metastable minimum, x+, and in 
the true minimum, X- , respectively. Eq. i|15|) was studied 
by Huse and Fisher in Ref. 0, where they showed that 
the asymptotic value of the residual energy e res (r) = 



FIG. 2: Same as Fig. for the asymmetric potential in 
Eq. JTJ} with a + = 1.25, a_ = 0.75 (dashed line in (a,b), 
inset of (c)). Notice the different behavior of RT and IT, in 
the present case. 

Ap + (r) is given by: 

£res (t) ~ const (71-) B (11171-) B , (16) 

where 7 is a constant. So, apart from logarithmic cor- 
rections, the leading power-law behavior is of the form 
£res ~ t -&v/b ^ wuere ^he exponent is controlled by the 
ratio Ay / B between the energy splitting of the two min- 
ima Ay and the barrier B. As shown in Figs. [He) an d 
[SJc) (solid lines through solid circles), the asymptotic be- 
havior anticipated by Eq. (|16fl fits nicely our CA resid- 
ual energy data (solid circles), as long as the logarithmic 
corrections are accounted for in the fitting procedure^. 




FIG. 3: Instantaneous eigenvalues (a) and ground state wave- 
functions (b) of the Schrodinger problem Hip — Elp for differ- 
ent values of F, for the symmetric potential in Eq. I14H with 
a+ = a- = a = 1. 



FIG. 4: Same as Fig. [3] for the asymmetric potential in 
Eq. 1141 with a+ = 1.25, a_ = 0.75. Notice the clear Landau- 
Zener avoided crossing in (a), indicated by the arrow and 
magnified in the inset. 



Obviously, we can make the exponent as small as we wish 
by reducing the linear term coefficient 5, and hence the 
ratio Ay/B, leading to an exceedingly slow classical an- 
nealing. 

The behavior of the QA evolution is remarkably dif- 
ferent. Starting from Fig. we notice that IT and RT 
evolutions give very similar residual energies, definitely 
faster decaying than the CA data, while the correspond- 
ing final wavefunction only slowly narrows around the 
minimum of the potential. Notice also the asymptotic be- 
havior of the residual energy, t res {r) oc r^ 1 / 3 , indicated 
by the dashed line in Fig. [TJc) : this rather strange expo- 
nent is simply the appropriate one for the Schrodinger an- 
nealing with a linear schedule within an harmonic poten- 
tial (the lower minimum valley, see Sec. Ill A|l . The asym- 
metric potential results, shown in Fig. [3 are even more 
instructive. The initial wavefunction squared \tj)(x,t = 
0)| 2 corresponds to a quite small mass (a large To = 0.5), 
and is broad and delocalized over both minima (solid 
line). As we start annealing, and if the annealing time r 
is relatively short - that is, if r < r c , with a character- 
istic time r c which depends on which kind of annealing, 
RT or IT, we perform - the final wavefunction becomes 
mostly concentrated on the wrong minimum, roughly cor- 
responding to the ground state with a still relatively large 



Ti < To (see Fig.0J). The larger width of the wrong valley 
is crucial, giving a smaller quantum kinetic energy con- 
tribution, so that tunneling to the other (deeper) min- 
imum does not yet occur. By increasing r, there is a 
crossover: the system finally recognizes the presence of 
the other minimum, and effectively tunnels into it, with a 
residual energy that, once again, decays asymptotically 
as e r es(T~) oc t -1 / 3 (dashed line in Fig. Hfc)). There is 
a characteristic annealing time r c - different in the two 
Scrodinger cases, RT and IT - above which tunneling 
occurs, and this shows up as the clear crossover in the 
residual energy behavior of both IT and RT, shown in 
Fig. He). 

These findings can be quite easily rationalized by 
looking at the instantaneous (adiabatic) eigenvalues 
and eigenstates of the associated time-independent 
Schrodinger problem, which we show in Fig. 0{a,b). 
Looking at the instantaneous eigenvalues shown in Fig. 
0Ja) we note a clear avoided-crossing occurring at T — 
T lz ~ 0.038, corresponding to a resonance condition be- 
tween the states in the two different valleys of the po- 
tential. For r > Tlz the ground state wavefunction is 
predominantly concentrated in the wider but metastable 
valley, while for T < T^z it is mostly concentrated on the 
deeper and narrower global minimum valley. In the full 
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time-dependent RT evolution, transfer to the lower valley 
is a Landau-Zener problem22*2i: the characteristic time 
t c for the tunneling event is given by tlz = hoiTo/2irA 2 , 
where a is the relative slope of of the two crossing 
branches as a function of T, 2 A is the gap at the avoided- 
crossing point, and Tq is the initial value of the anneal- 
ing parameter. (For the case shown in Fig. 0] we have 
2A = 0.0062, a = 2.3, hence t lz ~ 18980, see right- 
most arrow in Fig.[5Jc).) The Landau-Zener probability 
of jumping, during the evolution, from the ground state 
onto the "wrong" (excited) state upon fast approaching 
of the avoided level crossing is P ex = e~ T t TLZ , so that 
adiabaticity applies only if the annealing is slow enough, 
t > tlz- The IT characteristic time is smaller, in the 
present case, than the RT one. We will comment fur- 
ther on this point in Sec. IIII Al In a nutshell, the reason 
for this is the following. After the system has jumped 
into the excited state, which occurs with a probabil- 
ity P ex — e~ T / TLZ , the residual IT evolution will filter 
out the excited state; this relaxation towards the ground 
state is controlled by the annealing rate as well as by 
the average gap seen during the residual evolution. Nu- 
merically, the characteristic time r c seen during the IT 
evolution is of the order of H/(2A), see leftmost arrow in 
Fig. 0c), rather than being proportional to 1/A 2 as tlz 
would imply. 

Obviously, instantaneous eigenvalues/eigenvectors can 
be studied for the Fokker-Planck equation as well; their 
properties, however, are remarkably different from the 
Landau-Zener scenario just described for the Schrodinger 
case. Fig. |Sfc) shows the first four low-lying eigenvalues 
of the FP equation as a function of T (for a symmetric 
choice of the potential) , while Fig. [SJa,b) show the cor- 
responding eigenstates for two value of the temperature, 
T/V = 1 and T/Vq = 0.1. (The asymmetric potential 
cases are virtually identical, and are not shown). The 
lowest eigenvalue of the FP operator is identically and 
the corresponding eigenvector— is the Boltzmann distri- 
bution e~ v ( x ^ kBT , with roughly symmetric maxima on 
the two valleys. The first excited state correspond to 
distribution peaked on the two valley but with a node at 
the origin, and is separated from the ground state by an 
exponential small Arrhenius-like gap e ~ B / kBT . Higher 
excited states are separated by a very large gap, so that, 
effectively, only the two lowest lying states dominate the 
dynamics at small temperature. The reduction of a con- 
tinuum double- well FP classical dynamics onto a discrete 
effectively quantum two-level system, previously noticed, 
is quite evident from this form of the spectrum. On 
the contrary, the true quantum case never allowed for 
a discrete two-level system description whatsoever, ex- 
cept perhaps for large T. For small enough T < Tlz, in 
particular, the tower of oscillator states within the valley 
at x- is always very close in energy to the actual ground 
state, and the quantum annealing evolution reduces ef- 
fectively to a particle in a single harmonic well. This 
explains the rather large width of the final distributions 
P(x, t) observed in the quantum case. 
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FIG. 5: Instantaneous Fokker-Planck eigenvalues (panel (c)) 
as a function of temperature T, and the corresponding eigen- 
states for two values of T (panels (a) and (b)). The potential 
is here the symmetric one, V sym in Eq. ll3l with Vo = 1, a = 1, 
S = 0.1. Similar results (not shown) are obtained for the 
asymmetric double well potential Vkaym. 



Summarizing, we have found that QA and CA pro- 
ceed in a remarkably different way. CA is sensitive to the 
height of the barrier, more precisely to the ratio Ay jB 
between the energy offset Ay of the two minima, and 
the barrier height B. On the contrary, QA crucially de- 
pends on the tunneling probability between the two val- 
leys, which is reflected in a Landau-Zener (avoided cross- 
ing) gap: a wide tunneling barrier is obviously bad for 
QA. Finally, we noticed that RT and IT proceed with 
somewhat different characteristic times: we discuss this 
issue a bit more in the following section. 
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Real- versus imaginary-time Schrodinger 
evolution 



A Schrodinger dynamics in imaginary-time (IT) is 
clearly much more convenient than that in real-time (RT) 
for simulations on current classical computers, but it 
makes a difference in the final results? The answer to 
this question is, we believe: no, it does not make a dif- 
ference, in essence, although IT does give quantitatively 
better results. 

To qualify this statement, let us denote by ^^'(t)) 
the solution of the Scrodinger equation 



£^l* (e) (*)> = Wei 



H kin (t)]\^(t)) 



(17) 



where we assume that l^o) 1S the ground state of the 
initial Hamiltonian at time t$ , H c i + Hki n (to ) i while £ = 
ih for RT or £ = —IK for IT. By definition, the final 
residual energy after annealing up to time tf = to + r, 
where the kinetic energy is finally turned off, is given by: 



(r) 



(*M(t + T)\H cl \yM(t + T)) 

(*«)(t +r)|*(«(to+r)) 



E, 



Opt 



(18) 



We conjecture that the residual energies for the two al- 
ternative way of doing a Schrodinger evolution verify the 
following: i) the IT residual energy is not larger then the 
RT one, that is 



t {IT) (r) <^ rt Ht) , 

u res V / — res \ J ' 



(19) 



and ii) in many problems, the leading asymptotic be- 
havior, for t — * oo, might be identical for eres (t) and 

Expectation (i) seems very reasonable: it is simply 
inspired by the time-independent case, where it is well 
known that the IT Schrodinger dynamics tends to "filter 
the ground state" out of the initial trial wave function, as 
long as the gap between the GS and the first excited state 
is non-zero. However, we have here a time-dependent sit- 
uation, and the result is a priori not guaranteed. We do 
not have a proof of this statement, but we have veri- 
fied it in all the cases where an explicit integration of the 
Schrodinger equation has been possible (see, for instance, 
the results of the previous Section). (Needless to say, we 
have no proof of (ii) either, but, again, it never failed in 
all our tests.) 

The simplest time-dependent problem where one can 
test our conjectures, is the discrete two-level system 
(TLS) problem. Here, in terms of Pauli matrices, H c i = 
Aa z , while H kln (t) = -T{t)a x , with T(t) = -vt. The 
full H(t) is therefore 



H(t) = Aa z - T{t)a x 



(20) 



The annealing interpretation is very simple: the classical 



P (0) 




y = A /4v 



FIG. 6: The probability P ex (0) of ending up into the ex- 
cited state, given by Eg. 1211 for the discrete two-level system 
problem in Eq. 1201 for both imaginary-time (IT, dashed 
line) and real-time (RT, solid line) Schrodinger annealing. 
The large-7 behavior of P ex (0) is, in both cases, given by 
Pex(0) « l/(256 7 2 ). 



from the excited state | t) by a gap 2A. The kinetic 
term induces transitions between the two classical states. 
Starting from the ground state of H(to) at time to = — T 



we let the system evolve up to to time t 



f 



t + 



t = 0, 

when the Hamiltonian is entirely classical, H(tf = 0) = 
H c i = A<7 Z . The probability of missing the instantaneous 
final ground state | J.), ending up with the excited state 
| T), is: P«(0) = |(T |*^(0))| 2 /<* (O (0)|*«)(0))._ In 
principle, P ex depends, for given A, both on the initial 
r(io) = Tq = vt and on the annealing time r. The 
really important parameter, however, turns out to be the 
ratio v between these two quantities, which determines 
the "velocity of annealing": Taking t — > oo (i.e., to - ¥ 
— oo), and rh — > oo with T(t) = —vt for every t, the 
problem can be solved analytically (in terms of parabolic 
cylinder functions, see for instance Ref.0for the RT case) 
for both RT and IT. The probability P ex {fy of ending into 
the excited state can be expressed in terms of the variable 
7 = A 2 /4w. The explicit expressions, in terms of Gamma 
functions, are: 



R 



LR + ll 



2(l + |i?| 2 ) 

j_jXi+^)_ 
v^r(i/2 + ^) ' 



(21) 



optimal state is | J.), with energy E, 



opt 



-A, separated 



where 4>q — 3tt/4 and z = i^j for RT, while 4>o = ft and 
z = 7 for IT. A plot of P ex for both RT and IT is shown 
in Fig. as a function of 7 = A 2 /4t>. Note that: i) 
the IT-result for P ex (dashed line) is always below the 
RT-result, ii) the difference between the two curves is 
only quantitative: one can verify analytically that the 
leading behavior for large 7 is the same in both cases, 
i.e., P ex w l/(2567 2 ). Similar results are obtained by 
direct numerical integration of the Scrodinger equation 
for finite Tq and r, and with other forms of T(t). 
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FIG. 7: Comparison between the RT (solid lines) and the IT 
(dashed lines) evolution of a Landau-Zener problem, Eq. 1221 
for several values of the tunneling gap 2A (the values of 
A shown are A = 0.4, 0.2, 0.1, 10" 2 , 10" 3 , 10" 4 , 1(T 5 , 1(T 6 , 
while v — 1). The inset shows the two instantaneous eigen- 
values of the problem, E±(t), as a function of t. 



trivial effect of filtering, if on one hand explains the dis- 
crepancy between the IT and the RT evolution observed 




x/a 



FIG. 8: Parabolic washboard potential resulting in a loga- 
rithmically slow classical annealing. The minima are regularly 
located at positions Xi = ia, and the dashed line shows the 
parabolic envelope potential. 



With the same toy model, we can illustrate another 
point raised in the previous Section: what happens to the 
IT evolution after a Landau-Zener avoided crossing gap is 
encountered. The Hamiltonian we consider is essentially 
that in Eq. 1201 simply rotated in spin space, 

H(t) = -vto z - Aa x . (22) 

In absence of the tunneling amplitude A, the two en- 
ergy levels would cross at t = 0, while for A > 
the two instantaneous eigenvalues are simply E±(t) — 
±y/(vt) 2 + A 2 (see inset in Fig. [7J. Starting with the 
system in the ground state at t = — oo, we can monitor 
the probability of getting onto the excited states at any 
time t, which we plot in Fig. for both the RT and the 
IT evolution and for several values of A (taking v = 1). 
The RT data provide an illustration of the well-known 
Landau-Zener result: after a (relatively short) tunnel- 
ing time, and possibly a few oscillations, the probability 
of of getting onto the excited state saturates to a value 
given by P ex (t = oo) = e -^^l^. As for the IT data, 
the initial (tunneling) part and the subsequent plateau 
of the curves are similar to the RT case: the plateau 
value attained, call it P* x , is indeed very close to the 
RT saturation value (in fact, asymptotically the same 
for A — > 0); after that, the IT evolutions starts to filter 
out the ground state component - initially present in the 
state with a small amplitude 1 — P* x - through the usual 
mechanism of suppression of excited states, leading to a 
Pex{t) which is nicely fit by the curve 

P* e - a fo dt'[E+(t')-E^(t')] 

Pex{t) = (i _ p*) + p* e - 2 Io dt>\E + {t>)-E_{t>)] ' (23) 
which asymptotically goes to zero as t — > oo. This rather 



in the asymmetric double well case of the previous sec- 
tion, is, on the other hand, of no harm at all: on the 
contrary, it provides a quantitative improvement of IT 
over RT. 

In summary, the essential equivalence of IT and RT 
Schrodinger annealing (with, moreover, a quantitative 
improvement of IT over RT) justifies practical implemen- 
tations of quantum annealing based on imaginary-time 
Quantum Monte Carlo schemes. 

IV. ONE-DIMENSIONAL CURVED 
WASHBOARD: A POTENTIAL WITH MANY 
MINIMA 

After discussing at length the annealing problem in a 
potential with one minimum and with two minima, we 
wish to move on to a multi-minima problem, however 
simple. There are simple but interesting one-dimensional 
potentials which allow us to do that. The first example 
was proposed and solved by Shinomoto and Kabashima 
in Ref. and consists in a parabolically shaped wash- 
board potential. This example will display a logarithmi- 
cally slow classical annealing, showing CA may run into 
trouble even in simple models with no complexity what- 
soever, whereas quantum mechanics can do much better 
in this case. Consider a wiggly one-dimensional potential 
with barriers of individual height sa B separating differ- 
ent local minima, regularly located a distance a apart one 
from each other, i.e., at positions Xi = ai. The ith-local 
minimum is at energy Ci = ka 2 i 2 /2, so that the resulting 
envelope is parabolic. In order to study the dynamics of 
a particle in this potential, a good starting point is to 
write the master equation for the probability Pi(t) that 
the particle is in the ith-valley at time t: 
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-im) =P i+l {t)e~P B ^ +P l - 1 {t)e-^' -Pi(t) [ e -^.^+ c -^.«-] , (24) 
7 at 1 J 

I 



where 7 is an attempt frequency, Bij is the effective bar- 
rier from i to j, and /3 = l/kgT. This is a well justified 
starting point, in view of the results of the previous sec- 
tions fSecs. ITTTI and ITTA*) ! . showing that classical anneal- 
ing is extremely fast in reducing the width of a prob- 
ability distribution within each valley to an essentially 
delta-function like sequence of peaks of strength Pj . The 
actual form of the effective barriers depends on the way 
we model the details of the potential, with the only con- 
straint that detailed balance is satisfied, i.e., 

~ = £ i+i — e i — , 

in such a way that the stationary solution, for constant T, 
is simply the Boltzmann probability distribution, Pi (t — > 
00) oc exp(—ei/kBT). Ref. Il6l takes Sj+i^ = -Bj,i-i = B, 
while Bi_i + i = B + Ai and Sj-i^ = B + with 
Aj = e i+ i — 6j. The potential energy of the valleys enters 
only through the Sjj, which control the probability of 
making transitions between valleys. 

In order to study Eq. (|24|l , Shinomoto and Kabashima 
introduced a continuum limit, by defining a macroscopic 
coordinate x, such that the mimima are at Xi = ia, and 
writing the equation governing the probability P{x, t) in 
the limit a — > 0. The derivation involves writing Pi±\{t) 
in terms of derivatives of P(x,t), keeping consistently 
terms up to order a 2 and expanding exponentials with 
the assumption that ksT / {ka 2 ) » 1. The continuum 
limit equation governing the evolution of P{x, t) turns 
out to be a Fokker-Planck (FP) equation, Eq. (J3J, with 
an effective diffusion constant of the form 

D cS {T)= ia 2 e- B l k * T , (25) 

?/(T) = ksT/ D e g(T), and an effective drift potential 
V(x) = kx 2 /2 given by the macroscopic parabolic enve- 
lope potential. In order to study the annealing properties 
of the system, one can then follow exactly the same steps 
leading to Eq. J7J, which applies here too, except that 
now D t is substituted by D c ff(T(t)), which has an expo- 
nential activated behavior, D c g(T) oc e _s / fcBT . This 
exponentially activated D c g(T) changes the annealing 
behavior in a drastic way. Recall that the CA expo- 
nent Qca of Sec. IjH A|l decreases towards zero as the 
exponent an in the relationship D(T) oc T a ° increases. 
Since, close to T = 0, e ~ B/kBT « T a ° for any arbi- 
trarily large «c, we could suspect that the behavior of 
£res(T~) will no longer be a power law. In fact, the surpris- 
ing result of this exercise-^ is that the optimal annealing 
schedule T(t) is logarithmic and that e res (t) converges to 
at best as^ 

e res {t) ~ log(t)- 1 . (26) 



The reason behind this result is that the time derivative 
of e res (t) becomes exponentially small as one anneals T 
towards 0, due to the presence of e _s /' CsT in the diffusion 
constant, and an exponentially small derivative brings - 
not surprisingly - a logarithmically slow decrease of the 
function. To put it more physically, consider solving Eq. ( 
01 for a time independent T; the solution is trivially 

^res(^) = ^ relax _j ^xelax = 7! j 9^ ^ B 

2 2^ka l 

We observe that the solution converges to the equilib- 
rium (equipartition) value ksT/2 exponentially with a 
characteristic time i re iax which itself increases exponen- 
tially fast with decreasing T. As a result, the system 
will never be able to follow the decreasing T till the end 
of the annealing, by maintaining roughly the equilibrium 
value e P ot = ksT/2. Indeed, if we assume for instance 
T(t) = Tq(1 — t/r), the relaxation of the systems will 
cease to be effective - i.e., the system will fall out of 
equilibrium - at a time t* , and temperature T* — T(t*), 
at which i rc i ax ~ t, i.e., when ksT* B/logjr. The 
residual energy at this point cannot be smaller then the 
equipartition value ksT*/2, hence e res rs B/logjr as 
well. This freezing and falling out of equilibrium for clas- 
sical systems with barriers seems to provide an ubiqui- 
tous source of logarithms in classical annealing^. 

How would one tackle this annealing problem quan- 
tum mechanically? As the quantum analog of the master 
equation Eq. (|24|l . we propose studying the Schrodinger 
evolution governed by a tight-binding Hamiltonian with 
on-site energies £j and hopping matrix-elements between 
adjacent sites /ij,t+i 

H = ^ e i4 C i - ( C i Cl + 1 + C !+l C i) ' ( 27 ) 

i i 

The justification and possible limitations of this starting 
point, over the actual original continuum problem will be 
discussed at the end of the section. Here it is enough to 
consider that the hopping matrix-elements hij+i depend 
on tunneling through the barrier separating i from i + 1. 
The precise form of /ij^+i is likely to be inessential, in- 
cluding its energy (and hence site) dependence, for which 
we will assume the semi-classical (WKB) form: 

hi, i+ i~h = h (^y\-VW , (28) 

r = h 2 j (2ma 2 ) being simply related to the quantum con- 
finement energy of a particle of mass m in a valley of 
size ~ a. Here ho and Vh are energy parameters related 
to the details of the potential and of the barrier, which 
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will play little or no role. If the mass of the particle m 
(and hence T and h i i+ i) is kept constant the particle 
will explore the potential due to the kinetic term in the 
Hamiltonian: the correspondence between the quantum 
and the classical formulation is that F plays the role of 
T, hi t i+x plays the role of the classical transition prob- 
abilities, the ground state wavefunction |>I^ r:G ' S ' > | 2 at a 
given value of T (or, equivalently, of the hopping term 
h) plays the role of the classical equilibrium Boltzmann 
distribution p^ T,e9 - ) . The question, once again, is how to 
anneal T, by reducing it as a function of time, T(t), in 
such a way as to squeeze the wavefunction of the system 
so that the average potential energy 



Cpot(t) 



Eil*i(*)l a 



(29) 



is minimal. 

As it turns out, the continuum limit is once again use- 
ful. One goes to the continuum exactly as in the FP 
case by using x% = ai, writing ijj(xi,£) = ^i{i)/^/a and 
expanding everything to order a 2 . When written in first 
quantized form, the Hamiltonian for the quantum parti- 
cle in the macroscopic continuum coordinate x is simply 



H(t) = -T cS (t)V 2 + V(x) 



(30) 



where the coefficient of the Laplacian r o ff is the quantum 
counterpart of the classical effective diffusion constant 



T eS (t) = a 2 h(t) = a 2 h 



Vh 

r(t) 



1/4 



-Vv h /r(t) 



(31) 



and V(x) = kx 2 /2, as in the classical case. The con- 
tinuum limit quantum problem has therefore exactly the 
form we have considered in Sec. (jll A|l . except that T(t) 
in Eq. (|10[1 is now substituted by an effective Laplacian 
coefficient T c ff which has a highly non-linear, in fact ex- 
ponential, dependence on the annealing parameter T(t). 
We know, however, from Sec. (|II Ajl . that a non- linear 
behavior of the type (1 — t/r) ar for the Laplacian co- 
efficient leads to a power- law decrease of e res (r) with 
an exponent SIqa which is, remarkably, an increasing 
function of ar, approaching 1 as «r -* oo. Therefore, 
contrary to the classical case, where an exponentially ac- 
tivated behavior of the diffusion constant D c g is strongly 
detrimental to the annealing (turning a power law into 
a logarithm), here the exponential WKB-like behavior 
of r e g will do no harm at all. Indeed, we numerically 
integrated the relevant equation for B t , Eq. (|10|l with 
r e g in place of T, using the exponential WKB expres- 
sion Eq. |^nj for r c ff while annealing T(t) with a lin- 
ear schedule, T(t) = Ta(l — t/r). The integration was 
performed, as usual, with a fourth-order Runge-Kutta 
method, and was carried on up to time t = t, when the 
kinetic term in the Hamiltonian ceases to exist. The 
numerical results (not shown) have a clear power-law 
behavior for the quantum annealed (QA) final residual 



energy e res (t = r) ~ T _!iflA , with a power law expo- 
nent SIqa which is compatible with 1. Once again, the 
exponent appears to be insensitive to the choice of the 
type of quantum evolution (RT versus IT), although the 
numerical values of residual energies always respect the 
inequality e^(r) > e£T s ("r). 

Before ending this section, we would like to discuss 
briefly the reason for treating by tight-binding, Eq. 1271 
what was originally a continuum problem with a well 
defined potential landscape. As we learned from the 
double-well case, there is never a clear-cut discrete model 
(a discrete two-level system, in that case) describing in 
a complete way the continuum Schrodinger problem, in 
all stages of the annealing. Obviously, when the mass of 
the particle is very small, the tight-binding approxima- 
tion contained in Eq. 1271 is not particularly good, since 
more than one state per valley is generally important 
to describe the wavefunction accurately. As the mass 
of the particle increases, however, the tight-binding ap- 
proach gets more and more appropriate, until a further 
limitation appears: when the mass if very large, it is 
not legitimate to neglect excited states within, say, the 
central valley compared to the lowest states localized in 
metastable valleys. We can imagine that the ultimate 
behavior of e res (r), in the quantum case, will be actu- 
ally dominated by the rather trivial problem of squeez- 
ing the wavefunction in the lowest central minimum, with 
its characteristic power-law exponent (1/3, for instance, 
for a linear schedule). There is, however, an interme- 
diate region, between the very short r scale, where the 
full details of the potential are important, and the very 
long t scale, where the trivial squeezing mentioned above 
sets in, and where the tight-binding approximation rea- 
sonably predicts a power-law exponent (of order 1) for 

We believe that one of the important points that makes 
QA so different from CA in the present case is that the 
spectrum of the instantaneous eigenvalues of the quan- 
tum problem does not show any dangerous Landau-Zener 
avoided-crossing, and, correspondingly, the ground state 
wavefunction is always more peaked in the central valley 
than elsewhere. As in the two-level disorder in the 

width of the different valleys would change this result. 



V. ROLE OF DISORDER 

Despite their disarming simplicity, the three case stud- 
ies above turn out to be extremely informative in quali- 
fying the profound difference of QA from CA, and their 
surprising consequences. We expect that these results 
will be very important in understanding more realistic 
QA problems. Of course, the cases studied, although in- 
structive, do not possess the real ingredient which makes 
annealing difficult, both in CA and QA, i.e., some form 
of disorder in the distribution of the minima. We be- 
lieve, for instance, that even an irregular landscape with 
many minima, as the double-cosine potential V(x) — 



12 




-30 -20 -10 10 20 



FIG. 9: Double cosine potential V(x) = cos(2irx) + cos((l + 
■\ZE)ttx), showing an irregular landscape with many minima. 

V\ cos (2irx) + Vi cos (2nrx) (with r an irrational num- 
ber) shown in Fig. ^ would already change drastically 
the behavior of QA from a power-law to a logarithm. On 
quite general grounds, Anderson's localization^ would 
predict that wavefunctions are localized for a genuinely 
disordered potential for large enough mass (i.e., small 
enough kinetic energy bandwidth) in any D > 2 (this 
localization occurs for all value of the mass in D = 1,2). 
Therefore, quantum annealing should always, via a cas- 
cade of Landau-Zener events, end up into some localized 
state which has, a priori, nothing to do with our search 
of the actual potential minimum. 

A very simply illustration of the crucial role of disorder 
is given by the D = 1 disordered Ising ferromagnet: 

ff = -E J ^K- r E< (32) 

i i 

where Ji > are positive random variables in the inter- 
val [0, 1], and r is the transverse field inducing quantum 
fluctuations. Obviously, the ground state is the ferromag- 
netic state with all spins aligned up (or down). However, 
arbitrarily weak values of the Ji can pin domain walls 
between up and down ferromagnetic regions, with a very 
small energy cost 2J^. For a finite system with peri- 
odic boundary conditions domain walls appear in pairs, 
and separate sections of the system with alternating ] 
and { ferromagnetic ground states. Given two domain 
walls pinned at weak-J^ points a distance L> 1 apart, 
healing the system via single spin flip moves requires flip- 
ping L spins, which can be a formidable barrier to tun- 
nel through. The system will have a very slow anneal- 
ing (quantum, as well as classical) while showing, at the 
same time, no complexity whatsoever: simple disorder is 
enough. 

VI. SUMMARY AND CONCLUSIONS 

Summarizing, we have compared Schrodinger versus 
Fokker-Planck annealing in various simple cases of one- 



dimensional potential, in particular a double well poten- 
tial and a parabolic washboard potential with many min- 
ima but no disorder. In all cases the two annealings, 
quantum and classical, were seen to behave in a remark- 
ably different way. Classical annealing is influenced only 
by the height of the barriers surrounding the relevant 
minima, via Arrhenius-like thermal promotion over the 
barriers, with probability distributions which are quite 
localized on those minima. Quantum annealing is influ- 
enced by the structure of the eigenvalue spectrum of the 
problem - very small Landau-Zener tunneling gaps asso- 
ciated to large barriers are highly detrimental to it. In 
some cases, for instance in the case of the parabolic wash- 
board potential, quantum annealing can be much more 
effective then classical annealing, but, generally speak- 
ing, both strategies suffer whenever the potential land- 
scape is disordered. As an outcome of the discussion, it is 
quite clear that quantum annealing, although potentially 
useful and sometimes more convenient than classical an- 
nealing, is not capable, in general, of finding solutions of 
NP-complete problems in polynomial time: indeed, quite 
interestingly, even trivial optimization problems (trivial, 
obviously, only from the complexity point of view), like 
the one-dimensional ferromagnetic random Ising model, 
can show a very slow annealing behavior. 
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APPENDIX A: ANNEALING IN A PARABOLIC 
POTENTIAL 

1. Classical Annealing (Fokker-Planck) case. 

If is straightforward to find the solution of the Fokker- 
Planck Eq. when the the potential is harmonic, 
V(x) = kx 2 /2, and the initial condition is the Boltzmann 
distribution P(x,t = 0) oc exp (— fcx 2 /(2fcsTo)). Indeed, 
it is simple to verify that the Gaussian Ansatz 

P(x,t)^C t e~ BtX \ (Al) 

fulfills the initial condition (if B t= o = Bq = k/(2ksTo)) 
and solves the FP equation, as long as the two functions 
B t and Ct satisfy the following ordinary differential equa- 
tions: 

f B t = 2D t ( - 2B 2 ) 
\Ct/C t = D t ( 1 ^ m -2B t ) 
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The normalization constant Ct turns out to be irrelevant 
in calculating the average potential energy which we need 



tpot 



(t) 



/ dxV(x)P(x,t) 
J dxP{x,t) 



k 
*Bt 



T~d~ > (A3) 



and can be completely forgotten, since the equation for 
B t does not involve it. The equation for B t appears to be 
non-linear, but can immediately transformed into a linear 
equation by dividing up both sides by B 2 and recognizing 
that the correct variable to use is precisely 1/B t , or better 
yet, e po t(t). In terms of e po t(t) we can therefore write a 
linear equation of the form: 



d 

dl €pot{t) 



kD t l 



k B T(t) 



'-pot 



(*) 



(A4) 



the initial condition being simply given by the equipar- 
tition value e pot (t = 0) = (fc/4B ) = k B T /2. An 
alternative wayi& of deriving Eq. (| A4|) consists in tak- 
ing the derivative with respect to time of both sides of 
Eq. (|A3|) . using then the FP equation for dP/dt on the 
left hand-side of the ensuing equation, and finally inte- 
grating by parts the terms containing spatial derivatives 
of P (this procedure results in a closed differential equa- 
tion for e po t{t) only if the potential V(x) is harmonic). As 
every one-dimensional linear first-order differential equa- 
tion, Eq. (|A4|) can be integrated by quadrature, the so- 
lution being: 



tpot 



(*) 



tpot 



(t = 0) e~ {k/2kB) d v D v/ T (y) 



+k / d^ e -(V3*a)#«WT(v).(A5) 



If we now anneal the temperature down to in a time r 
in the usual way, T(t) = Tq(1 — t/T) aT , assuming that the 
diffusion constant behaves as D(T) = D (T/T ) aD , we 
readily get an analytic expression for the residual energy 
at the end of the annealing, e res (T) = e po t(t = r), which 
can be shown to behave, for large r, as a power law: 



£re S (T) 



-Oca 



a 



CA 



ax(a D - 1) + 1 



(A6) 



Quite evidently, annealing proceeds here extremely fast, 
with a power-law exponent tie A that can increase with- 
out bounds (for instance if ap = 1) upon increasing the 
exponent olt of the annealing schedule. Notice, however, 
that large values of an are, on the contrary, detrimental. 



2. Quantum Annealing (Schrodinger) case. 

Consider now the problem of a particle moving in 
a parabolic potential V(x) = kx 2 /2, with a time- 
dependent mass, such that the Hamiltonian is given by: 



H{t) 



r(*)v 2 , 



(A7) 



where T(t) — h 2 /2m(t) denotes the coefficient of the 
Laplacian operator. The Schrodinger evolution of the 
wavefunction ^>(a;,i) is then, 



£dtil){x,t) = H(t)ip(x,t) , 



(A8) 



where £ = ih for a real time (RT) evolution, and £ = — h 
for an imaginary time (IT) evolution. Now, whereas gen- 
eral solutions of the time-dependent Schrodinger equa- 
tion for arbitrary initial condition ij){x,t = 0) are not 
easy, it turns out that if V(x) is quadratic, then any ini- 
tial Gaussian wavefunction propagates into a Gaussian, 
which is enough for our goal. In detail, write the follow- 
ing Ansatz for tp(x, t): 

ip(x, t) = C t e- BtX ^ 2 Real(Bt) = $t(B t ) > . 

(A9) 

Substituting the Ansatz for ij]{x,t) into the Schrodinger 
evolution (in RT or in IT), one immediately verifies that 
ip(x,t) satisfies Eq. (|A8() for arbitrary T(t) as long as B t 
and Ct satisfy the following ordinary differential equa- 
tions: 



-£B t =k- 2T(t)B 2 
<C t /C t = -T(t)B t 



(A10) 



Once again, the normalization constant C't turns out to 
be irrelevant in calculating the average potential energy 



£pot(t) 



JdxV(x)\ip(x,t)\ 2 
fdx\i;(x,t)\2 



(All) 



and can be completely forgotten, since the equation for 
Bt does not involve it. The initial condition Bt=o is, 
by assumption, such that at t = the system is in the 
ground state corresponding to To = T(t = 0). Such a 
ground state value Bq is easily calculated by equating to 
zero the right-hand side of Eq. (|A10|) with T(t) — » r , 
i.e., B = y/k/(2T ). 

A few general considerations can be based on a purely 
qualitative analysis of Eq. i|A10|) . Consider first the IT 
case, where the equation for B t reads HB t = k — 2T(t)B 2 . 
One can easily get convinced that B t is forced to be a real, 
positive and monotonically increasing function of t, i.e., 
B t > and B t > for t > 0. Therefore, we can easily 
write an inequality of the form: 

hB t = k- 2T(t)B 2 < k , 

from which we immediately conclude, by integrating over 
t the two sides of the inequality, that 



B T 



kr 

B < -r 

n 



i.e., the residual energy e res (t) = k/(4B T ) cannot de- 
crease faster than 1/t for r — * oo. 

We will now assume, without great loss of generality, 
that the Laplacian coefficient T(t) is given by T(t) = 
ro/(t/r), where r is an annealing time-scale (for in- 
stance the annealing time, when a linear schedule is used) 
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and f(t') is a positive decreasing function for t' > 
such that f{t' < 0) = 1. It is useful to switch to di- 
mensionless variables by measuring times in unit of r, 
t' = t/r, and B t in units of its initial ground state value 
B t= a = Bq. The appropriate dimensionless quantity is 
therefore b(t'; r) = Bt/Bg, with t' — t/r, where the para- 
metric dependence on the annealing time-scale r has been 
explicitly indicated. The equation for b(t'; t) is given by: 



6(t';r) 
6(0; t) 



a[l 
1 . 



}{t')b 2 {t'-r)] 



T^2kT /(-0 
(A12) 



where the dot, from now on, will denote a derivative with 
respect to t' . Notice that the parametric dependence 
on r is all buried into the coefficient a, which reabsorbs 
also the — £ appearing in the dynamics (RT versus IT). 
This kind of non-linear differential equation is of the well 
known Riccati form. It can be transformed into a lin- 
ear second-order differential equation by operating the 
following substitution 



b(t';r) = 



Mr) 



af(t')y(t';r) ' 



(A13) 



where, evidently, y is defined up to an overall normaliza- 
tion constant. Indeed, simple algebra shows that we can 
re-express Eq. (|A12|I as a second-order linear equation 
for y, as follows: 

f(t')y(t'; t) - j{t')y(t'-T) = a 2 f 2 (t')y(t'; r)(A14) 
y(0; r) = ay(0;r) . 



As long as f(t') ^ 0, it is simple to verify that the second- 
order equation for y(t'; r) can be also equivalently written 



d ( y{t';T) , 2 , 



dt 



fit') 



(A15) 



Finally, denoting by Y(t';T) the indefinite integral of 
y(t';T), such that Y = y, and integrating over t' both 
sides of Eq. (|A15|) . we can also write: 



Y{t'-r)=a 2 f{t')Y{t';T) 



(A16) 



Eq. (|A16p is easily solved when the annealing schedule 
for T(f) is linear, T(t) = T (l - t/r), i.e., when /(*') = 
(1 — t') av with a>r = 1, a case in which Eq. (|A16|) is of the 
Airy type. In the latter case, it is useful to perform a final 
change of independent variable to z = £ 2 |a|3(l — t'), so 
that, defining F(z) — Y(t';r), we can write the equation 
for F in the standard Airy form: 



(A17) 



The general solution of Eq. (|A17|I is given in terms of the 
two Airy's functions Ai(z) and Bi(z), 



F(z) = pAi(z) + jBi(z) 



(A18) 



where (3 and 7 are two constant coefficients. Going back 
to Y(t';r) and y(t';T), we then have the explicit expres- 
sions: 



Y(t';r) = /3^ 2 |«|i(l-i'))+7^(e 2 |«l i (l-0) 

Y(t'-r) = y(t';T) = -pe\a^At'(e\a\Hl~t'))~ 1 e\^B i '(e\a\ 2 Hl-t')) 

Y(t';r) = y{t'-T)=p\a^M'\e\^{l-t'))+ 1 \a\iBi"{e\a\i{l-t')) 

= a 2 (l-t'){/3^(e 2 Hi(l-t'))+7^(e 2 |«l l (l-0)} , (A19) 



where the prime indicates a derivative with respect to z, 
and we have used the property of Airy's functions that 
Ai"(z) — zAi(z) and Bi"(z) = zBi(z). Finally, substi- 



tuting back the expressions in Eq. I|A19|I into the original 
function B t , see Ea. lA13l we get: 



Bt — Br 



yith 



ar(l - t/T)y(t/T\ 



tiB \a\ 



(3Ai(e\a\i(l - t/r)) + yBi(f\a\i(l - t/r)) 
/W(£ 2 |a|i(l - t/T))+ 1 Bi'{e\a\i{\ - t/r)) 



(A20) 



r 



This general solution correctly depends on one parameter generality. We impose the initial condition B t= o = Bq 
only, i.e., 7//3, so that we can put (3 = 1 without loss of 
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by requiring: 

.1 A l {e\a\^)+ 1 Bi(e\^] 



= i 



Ai'(Z 2 \a\i) + jBi'(£ 2 \a\i] 
Solving for 7, we get: 

£ 2 H^(£ 2 |a|§)-,4z'(£ 2 |a|§; 



7 



^\a\iBi(^\a\i) - Bi'(£, 2 \a\i) 



(A21) 



(A22) 



Due to the asymptotic properties of the Airy's functions 
(Ai(z) — > 0, Ai'(z) — > and Bi(z) — > 00, B«'(z) — > 00, 
when 2: — > 00 on the real axis), we conclude that for the 
IT case: 



7^0 
so that, finally, 

B(r) = 



for a (or r) — > 00 



(A23) 



-Sola 



Ai(0)) +7-Bi(0) 



Ai'(0) 

/2fcr r N 



h75i'(0) 

r(i) 



s*r(|) 



(A24) 



where we used that Ai(0) = 3" 2 / 3 /r(2/3) and Ai'(0) = 
-3" 1 / 3 /r(l/3). On the other hand, for the RT case we 
have to take the limit z — ► — 00 (on the real axis) in- 
stead, and in that region all the Airy's functions oscil- 
late. Nevertheless is possible to show that the value of 7 
is uniformly bounded for r — > 00. 

We conclude, therefore, that for large r and with a 
linear annealing schedule, otr = 1, 



B(t) oc t~< 



(A25) 



and, consequently, the residual energy behaves asymp- 
totically as 



^res \T) 



4St{B T ) 



OC T 3 



(A26) 



The generalization of this result to an arbitrary annealing 
exponent ar > in T(t) = Tq(1 — t/r) ar , is a bit more 
involved. It is however possible to establish a generaliza- 
tion of Eq. (|A26|) . for arbitrary ar, that we checked by 
means of direct numerical integration. It reads: 



^res (^~) 



k 



■ 



4£(B T )~' QA « r + 2' 

(A27) 

an expression that holds true for both RT and IT anneal- 
ing. 



APPENDIX B: CLASSICAL ANNEALING WITH 
QUANTUM TOOLS: IMAGINARY TIME 
SCHRODINGER EVOLUTION OF THE 
FOKKER-PLANCK EQUATION. 

A side issue, but nonetheless an interesting one we wish 
to discuss here relates to classical annealing, and con- 



cerns the well-known mapping of a Fokker-Planck equa- 
tion onto an imaginary-time Schrodinger problem 2 ^, and 
its implications on the relationship between CA and QA. 
The bottom-line will be that the mapping does not im- 
ply that a FP-based CA is actually equivalent to QA, and 
moreover that such a mapping is not particularly useful 
in our annealing context. 

Consider, once again, the FP problem with a time- 
dependent temperature T(t) 



0_ 

dt 



P(x,t) 



1 



div (PVV) + AV 2 P 



(Bl) 



where both the friction coefficient and the diffusion con- 
stant are now time-dependent quantities, which we in- 
dicate by rjt = r}(T(t)) and D t = D{T(t)). In order 
to map the problem in Eq. IB1I onto an imaginary-time 
Schrodinger problem, the standard procedure 2 * is to pose 
P(x,t) — <&o(x,t)i/}(x,t) and to determine <&o(x,t) in 
such a way as to eliminate the non-Schrodinger-looking 
drift term, turning it onto a standard potential term. 
The algebra is trivial. One can show that the drift term 
is eliminated if, and only if, the $0 satisfies the equation: 

V$oOM) = -5— =j-$o(a:,t) , 
2r] t D t 

whose solution is readily found to be: 

$o(x,t) =C(t)e~ v/ ^ tDt) , (B2) 

with C(t) a function of time only, which can even be 
taken to be constant without loss of generality. By plug- 
ging P = $oip in the FP equation IB1I with $0 as above, 
one can show that the resulting equation for ip(x,t) is 
indeed of the Schrodinger form 



-D t V 2 ip(x, t) + V FP (x, t)ip{x, t) , (B3) 



with an effective potential Vpp given by 



V F p(x,t) = ^- 



(W) 2 



V 2 V 



d t $o{x,t) 
<P (x,t) 



(B4) 



The first term in Vpp is the standard effective poten- 
tial of the Riccati form obtained in the time-independent 
case24. The second piece in Vpp is absent in the time- 
independent case, and can be easily seen to be 2 ^: 

d t ® {x,t) = d_ f 1 

$ (z,t) 1 'dt \2 Vt D t 

The main point we want to stress is that, by annealing 
T(t) and hence D t , we not only reduce the coefficient of 
the Laplacian in Eq. IB3I but we also strongly modify the 
potential Vfp, at variance with a genuine QA where only 
the kinetic term is annealed down. The modifications of 
the potential are so strong that, at low temperature, the 
instantaneous eigenvalue spectrum associated to the FP 
equation, as discussed in Sec. IIIII is vastly different from 
that of the quantum double well system. 
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